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The density-matrix renormalization group (DMRG) applied to transfer ma- 
trices allows it to calculate static as well as dynamical properties of one- 
dimensional quantum systems at finite temperature in the thermodynamic 
limit. To this end the quantum system is mapped onto a two-dimensional 
classical system by a Trotter-Suzuki decomposition. Here we discuss two dif- 
ferent mappings: The standard mapping onto a two-dimensional lattice with 
checkerboard structure as well as an alternative mapping introduced by two of 
us. For the classical system an appropriate quantum transfer matrix is defined 
which is then treated using a DMRG scheme. As applications, the calculation 
of thermodynamic properties for a spin- 1/2 Heisenberg chain in a staggered 
magnetic field and the calculation of boundary contributions for open spin 
chains are discussed. Finally, we show how to obtain real time dynamics from 
a classical system with complex Boltzmann weights and present results for 
the autocorrelation function of the XXZ-cha.m. 



1 Introduction 

Several years after the invention of the DMRG method to study ground- 
state properties of one-dimensional (ID) quantum systems [1], Nishino showed 
that the same method can also be applied to the transfer matrix of a two- 
dimensional (2D) classical system hence allowing to calculate its partition 
function at finite temperature [2] . The same idea can also be used to calculate 
the thermodynamic properties of a ID quantum system after mapping it to 
a 2D classical one with the help of a Trotter- Suzuki decomposition [3, 4, 5]. 
Bursill et. al. [6] then presented the first application but the density matrix 
chosen in this work to truncate the Hilbert space was not optimal so that the 
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true potential of this new numerical method was not immediately clear. This 
changed when Wang and Xiang [7] and Shibata [8] presented an improved 
algorithm and showed that the density-matrix renormalization group applied 
to transfer matrices (which we will denote as TMRG from hereon) is indeed 
a serious competitor to other numerical methods as for example Quantum- 
Monte-Carlo (QMC). Since then, the TMRG method has been successfully 
applied to a number of systems including various spin chains, the Kondo 
lattice model, the t — J chain and ladder and also spin-orbital models [9, 10, 
11, 12, 13, 14, 15, 16, 17] . 

The main advantage of the TMRG algorithm is that the thermodynamic 
limit can be performed exactly thus avoiding an extrapolation in system size. 
Furthermore, there are no statistical errors and results can be obtained with 
an accuracy comparable to T = DMRG calculations. Similar to the T = 
DMRG algorithms, the method is best suited for ID systems with short 
range interactions. These systems can, however, be either bosonic or fermionic 
because no negative sign problem as in QMC exists. Most important, there 
are two areas where TMRG seems to have an edge over any other mimerical 
methods known today. These are: (1) Impurity or boundary contributions, and 
(2) real-time dynamics at finite temperature. As first shown by Rommer and 
Eggert [18], the TMRG method allows it to separate an impurity or boundary 
contribution from the bulk part thus giving direct access to quantities which 
are of order 0{1/L) compared to the 0{1) bulk contribution (here L denotes 
the length of the system). We will discuss this in more detail in section 5. 
Calculating numerically the dynamical properties for large or even infinite ID 
quantum systems constitutes a particularly difficult problem because QMC 
and TMRG algorithms can usually only deal with imaginary-time correlation 
functions. The analytical continuation of numerical data is, however, an ill- 
posed problem putting severe constraints on the reliability of results obtained 
this way. Very recently, two of us have presented a modified TMRG algorithm 
which allows for the direct calculation of real-time correlations [19]. This new 
algorithm will be discussed in section 6. 

Before coming to these more recent developments we will discuss the def- 
inition of an appropriate quantum transfer matrix for the classical system in 
section 2 and describe how the DMRG algorithm is applied to this object in 
section 3. Here wc will follow in parts the article by Wang and Xiang in [20] 
but, at the same time, also discuss an alternative Trotter-Suzuki decomposi- 
tion [15, 16]. 

2 Quantum transfer matrix theory 

The TMRG method is based on a Trotter-Suzuki decomposition of the parti- 
tion function, mapping a ID quantum system to a 2D classical one [3, 4, 5]. In 

the following, we disc;uss both the standard mapping introduced by Suzuki [5] 
as well as an alternative one [15, 16] starting from an arbitrary Hamiltonian 
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i? of a ID quantum system with length L, periodic boundary conditions and 
nearest neighbor interaction 

L 

H = ^h,,,+i. (1) 

i=l 

The standard mapping, widely used in QMC and TMRG calculations, is de- 
scribed in detail in [20]. Therefore we only summarize it briefly here. First, 
the Hamiltonian is decomposed into two parts, H = He + Hg, where each part 
is a sum of commuting terms. Here (Ho) contains the interactions 
with i even (odd). By discretizing the imaginary time, the partition function 
becomes 

Z = Tre-f^" = \im TV { [e-^^^e-^^"]"^} (2) 

with e — P/M, l3 being the inverse temperature and M an integer (the so 
called Trotter number). By inserting 2M times a representation of the identity 
operator, the partition function is expressed by a product of local Boltzmann 
weights 

^k,k+i — vk^k 1 6 i-'fe+i*fe+i/ ' y-^) 

denoted in a graphical language by a shaded plaquette (see Fig. 1). The sub- 
scripts i and k represent the spin coordinates in the space and the Trotter 
(imaginary time) directions, respectively. A column-to-column transfer matrix 
7m, the so called quantum transfer matrix (QTM), can now be defined using 
these local Boltzmann weights 

7m = (ti,2T3,4 • • • T2M-1,2m) (t2,3T4,5 • • • T2M,l) • (4) 

and is shown in the left part of Fig. 1. The partition function is then simply 
given by 

Z = Tv T^'^ . (5) 

The disadvantage of this Trotter-Suzuki mapping to a 2D lattice with checker- 
board structure is that the QTM is two columns wide. This increases the 
amount of memory necessary to store it and also complicates the calculation 
of correlation functions. 

Alternatively, the partition function can also be expressed by [15, 16] 

Z= Jim Tr{[Ti(e)T2(e)]^/n , (6) 

with Ti,2(e) = Tr^l exp[-ei7 + ©(e^)]. Here, Tr,l are the right- and left-shift 
operators, respectively. The obtained classical lattice has alternating rows 
and additional points in a mathematical auxiliary space. Its main advantage 
is that it allows to formulate a QTM which is only one column wide (see right 
part of Fig. 1). The derivation of this QTM is completely analogous to the 
standard one, even the shaded plaquettes denote the same Boltzmann weight. 
Here, however, these weights are rotated by 45° clockwise and anti-clockwise 
in an alternating fashion from row to row. Using this transfer matrix, 7m, the 
partition function is given by Z = Tr 7^ . 
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Fig. 1. The left part shows the standard Trotter-Suzuki mapping of the ID quan- 
tum chain to a 2D classical model with checkerboard structure where the vertical 
direction corresponds to imaginary time. The QTM is two-column wide. The right 
part shows the alternative mapping. Here, the QTM is only one column wide. 



2.1 Physical properties in the thermodynamic Hmit 

The reason why this transfer matrix formalism is extremely useful for numer- 
ical calculations has to do with the eigenspectrum of the QTM. At infinite 
temperature it is easy to show [21] that the largest eigenvalue of the QTM 
7m (Tm) is given by S*^ (S) and all other eigenvalues are zero. Here S de- 
notes the number of degrees of freedom of the physical system per lattice 
site. Decreasing the temperature, the gap between the leading eigenvalue Aq 
and next-leading eigenvalues An {n > 0) of the transfer matrix shrinks. The 
ratio between Aq and each of the other eigenvalues yl„, however, defines a 
correlation length l/Cn = lii|^o/^n| [20, 21]. Because a ID quantum system 
cannot order at finite temperature, any correlation length ^„ will stay finite 
for T > 0, i.e., the gap between the leading and any next-leading eigenvalue 
stays finite. Therefore the calculation of the free energy in the thermodynamic 
limit boils down to the calculation of the largest eigenvalue Aq of the QTM 
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Here the interchangeability of the hmits L ^ oo and e ^ has been used 
[5] . Local expectation values and static two-point correlation functions can be 
calculated in a similar fashion (sec e.g. [20] and [21]). In the next section, we 
are going to show how the eigenvahies of the QTM are computed by means of 
the density matrix renormalization group. This is possible since the transfer 
matrices are built from local objects. Instead of sums of local objects we 
are dealing with produc;ts, but this is not essential to the numerical method. 
However, there are a few important differences in treating transfer matrices 
instead of Hamiltonians. At first sight, these differences look technical, but at 
closer inspection they reveal a physical core. 

The QTMs as introduced above are real valued, but not symmetric. This 
is not a serious drawback for numerical computations, but certainly inconve- 
nient. So the first question that arises is whether the transfer matrices can be 
symmetrized. Unfortunately, this is not the case. If the transfer matrix were 
replaceable by a real symmetric (or a hermitean) matrix all eigenvalues would 
be real and the ratios of next-leading eigenvalues to the leading eigenvalue 
would be real, positive or negative. Hence all correlation functions would 
show commensurability with the lattice. However, we know that a generic 
quantum system at sufficiently low temperatures yields incommensurate os- 
cillations with wave vectors being multiples of the Fermi vector taking rather 
arbitrary values. 

Therefore we know that the spectrum of a QTM must consist of real 
eigenvalues or of complex eigenvalues organized in complex conjugate pairs. 
This opens the possibility to understand the QTM as a normal matrix upon 
a suitable choice of the underlying scalar product. Unfortunately, the above 
introduced matrices are not normal with respect to standard scalar products, 
i.e. we do not have [Tm,T^] = 0. 



Next, we describe how to increase the length of the transfer matrix in imagi- 
nary time, i.e. the inverse temperature, by successive DMRG steps. Like in the 
ordinary DMRG, we first divide the QTM into two parts, the system S and 
the environment block E. Using the QTM, 7m, the density matrix is defined 



modynamic limit. As in the zero temperature DMRG algorithm, a reduced 
density matrix ps is obtained by taking a partial trace over the environment 



3 The Method - DMRG algorithm for the QTM 




ps = TrB{K)K|} • 



(9) 
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Note that this matrix is real but non-symmetric, which compHcates its numer- 
ical diagonalization. It also allows for complex conjugated pairs of eigenvalues 
which have to be treated separately (see [21] for details). 

In actual computations, the Trotter-Suzuki parameter e is fixed. Therefore 

the temperature T ~ 1/eM is decreased by an iterative algorithm M M+1. 
In the following, the blocks of the QTM, Tm, are shown in a 90°-rotated view. 

1. First we construct the initial system block F consisting of M plaquettes 
so that 5^ < iV < S^+S where S is the dimension of the local Hilbert 
space and N is the number of states which we want to keep. ns,n'g are 



ns 




n^ 

Fig. 2. The system block F. The plaquettes are connected by a summation over 
the adjacent corner spins. 

block-spin variables and contain N = states. The ■ TV^-dimensional 
array F{a, Ug, r, n'g) is stored. 
2. The enlarged system block F{a,ns,S2,T,s'2,n'g), a S"' • iV^-dimensional 
array, is formed by adding a plaquette to the system block. If hi^i+i is real 
and translationally invariant, the environment block can be constructed by 
a 180°-rotation and a following inversion of the system block. Otherwise 
the environment block has to be treated separately like the system block. 
Together both blocks form the superblock (see Fig. 3). 




Fig. 3. The superblock is closed periodically by a summation over all a states. 



3. The leading eigenvalue Aq and the corresponding left and right eigenstates 

(iZ^o^l =>f^(si,n„S2,ne) , \^^) =^'^{s[,n'^,4,n'J (10) 

are calculated and normalized {'l^if I'I'q ') = 1- Now thermodynamic quan- 
tities can be evaluated at the temperature F = l/(2e(M + 1)). 



DMRG applied to transfer matrices 7 

4. A reduced density matrix is calculated by performing the trace over the 
environment 

Ps{n',,4\n,,S2)=J2K)('^o\ (11) 

sl,ne 

= ^ iZ'^(si,n^,S2,ne)!f^(si,ns,S2,ne) 

sl,ne 

and the complete spectrum is computed. A A^x (5-iV)-matrix V^{ns\ns, S2) 
{V^{n'Jn'g,s'2)) is constructed using the left (right) eigenstates belonging 
to the N largest eigenvalues, where (n^) is a new renormalized block- 
spin variable with only A^ possible values. 

5. Using and the system block is renormalized. The renormalization 



Us 




lis 

Fig. 4. The renormalization step for the system block. 



is given by 

r{(7,ns,T,n'^) = ^ ^ V^{ns\ns,S2)r{cr,ns,S2,T,s'2,n'^)V^{n'Jn'^,S2) . 

(12) 

Now the algorithm is repeated starting with step 2 using the new system 
block. However, the block-spin variables can now take A^ instead of TV values. 

4 An example: The spin-1/2 Heisenberg chain with 
staggered and uniform magnetic fields 

As example, we show here results for the magnetization of a spin-1/2 Heisen- 
berg chain subject to a staggered magnetic field hg and a uniform field 
hu = qij-bH/J 

'H = JY,[^^^^+r-KSt-{-lyKS^] , (13) 
i 

where H is the external uniform magnetic field and g the Landc factor. An 
effective staggered magnetic field is realized in spin-chain compounds as for 
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example copper pyrimidine dinitrate (CuPM) or copper benzoate if an exter- 
nal uniform magnetic field H is applied [22] . For CuPM the magnetization as 
a function of applied magnetic field H has been measured experimentally. In 
Fig. 5 the excellent agreement between these experimental and TMRG data 
at a temperature T = 1.6 K with a magnetic field applied along the c" axis 
is shown. Along the c" axis the effect due to the induced staggered field is 
largest (see [23] for more details). Note that at low magnetic fields the TMRG 
data describe the experiment more accurately than the exact diagonalization 
(ED) data, because there are no finite size effects (see inset (a) of Fig. 5). For 




I^qH (T) 

Fig. 5. TMRG data (solid lino) and experimental magnetization curves (circles) 
for CuPM at a temperature T" = 1.6 K with the magnetic field applied along the 
c" axis. For comparison ED data for a system of 16 sites and T = are shown 
(dashed lines). Hero J/fc_B = 36.5 K, hu = gusH/J, hs — O.llhu and g — 2.19. 
Inset (a): Magnetization for small magnetic fields. Inset (b): Susceptibility as a 
function of temperature T at H = calculated by TMRG. 

a magnetic field H applied along the c" axis a gap, A oc iJ^/"^, is induced with 
multiplicative logarithmic corrections. For — > and low T the susceptibility 
diverges x ~ I/T because of the staggered part [24] (see inset (b) of Fig. 5). 
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5 Impurity and boundary contributions 

In recent years much interest has focused on the question how impurities 
and boundaries influence the physical properties of spin chains [25, 26, 27, 
28, 29]. The doping level p defines an average chain length L = 1/p — 1 
and impurity or boundary contributions are of order ~ 0{1/L) compared to 
the bulk. This makes it very difficult to separate these contributions from 
finite size corrections if numerical data for finite systems (e.g. from QMC 
calculations) are used. TMRG, on the other hand, allows to study directly 
impurities embedded into an infinite chain [18] . We will discuss here only the 
simplest case that a single bond or a single site is different from the rest. The 
partition function is then given by 



T^M '^impj ) (14) 

where Ti^p is the QTM describing the site impurity or the modified bond. In 
the thermodynamic limit the total free energy then becomes 

F = -TlnZ = LU^,^ + = -iTlnylo - T\n{\^JAo) , (15) 

with Aq being the largest eigenvalue of the QTM, 7m, and Ximp = ^Q\%mpWQ)- 

As example, we want to consider a semi-infinite spin-1/2 XXZ-oham with 
an open boundary. In this case translational invariancc is broken and field 
theory predicts Friedel-type oscillations in the local magnetization {S^(r)) 
and susceptibility %('') = d{S^{r))/dh near the boundary [30, 31]. Using the 
TMRG method the local magnetization can be calculated by 

(S-W) = i^^\nsj'-^r...\K) , ,,,, 

^^Q^imp 

where T{S^) is the transfer matrix with the operator included and is 

the transfer matrix corresponding to the bond with zero exchange coupling. 
Hence %^p\<F'^) is nothing but the state describing the open boundary at the 
right. In Fig. 6 the susceptibility profile as a function of the distance r from 
the boundary for various temperatures as obtained by TMRG calculations 
[31] is shown. For more details the reader is referred to [18] and [31]. 



6 Real time dynamics 

Finally, we want to discuss a very recent development in the TMRG method. 
The Trotter-Suzuki decomposition of a ID quantum system yields a 2D clas- 
sical model with one axis corresponding to imaginary time (inverse temper- 
ature) . It is therefore straightforward to calculate imaginary-time correlation 
functions (CFs). Although the results for the imaginary-time CFs obtained by 
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Fig. 6. Susceptibility profile for A — 0.6 and different temperatures T. N = 240 
states have been kept in the DMRG algorithm. The lines are a guide to the eye. 

TMRG are very accurate, the results for real-times (real-frequencies) involve 
errors of unknown magnitude because the analytical continuation poses an ill- 
conditioned problem. In practice, the maximum entropy method is the most 
efficient way to obtain spectral functions from TMRG data. The combination 
of TMRG and maximum entropy has been used to calculate spectral functions 
for the XXZ-chain [17] and the Kondo-lattice model [14]. However, it is in 
principle impossible to say how reliable these results are because of the afore 
mentioned problems connected with the analytical continuation of numerical 
data. It is therefore desirable to avoid this step completely and to calculate 
real-time correlation functions directly. 

A TMRG algorithm to do this has recently been proposed by two of us 
[19]. Starting point is an arbitrary two-point CF for an operator Or{t) at site 
r and time t 



Here we have used the cyclic invariance of the trace and have written the 
denominator in analogy to the numerator. In the following we will use the 
standard Trotter-Suzuki decomposition leading to a two-dimensional checker- 
board model. 




(17) 
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The crucial step in our approach to calculate real-time dynamics directly 
is to introduce a second Trotter- Suzuki decomposition of cxp{—iSH) with 
S = t/N in addition to the usual one for the partition function described in 
section 2. We can then define a column-to-column transfer matrix 

T2N,M — (ti,2T3,4 • • • T2M-1,2m)(t2,3'''4,5 ■ • • 7'2M,2M+l) (18) 
(w2M+l,2M+2 • • • V2M+2N-1,2M+2n){v2M+2,2M+3 ' ' ' V2M+2N ,2M+2N +l) 
{V2M+2N+\,2M+2N+1 • ' ■ V2M+iN-l,2M+iN){V2M+2N+2,2M+2N+S ' ' ■V2M+iN,l) 

where the local transfer matrices have matrix elements 

t(4«^'I4+i5^) = (44+V-^'-^M4+is^) (19) 

and V is the complex conjugate. Here i = is the lattice site, 

k = 1, ■ • ■ , 2M (/ = 1, • • • , 2N) the index of the imaginary time (real time) 
slices and s\g^^ denotes a local basis. The denominator in Eq. (17) can then 

L 12 

bo represented by Tr(7^y^^^) where N,M,L ^ oc. A similar path-integral 
representation holds for the numerator in (17). Here we have to introduce an 
additional modified transfer matrix T2n,m (O) which contains the operator O 
at the appropriate position. For r > 1 we find 

Tr(Tf6)T[''/21-iT(d)T^/2-I''/21-i) 
mtmO)) = ^ hm^^^lim ^ ^ ^ ^^^^12] ^ 

= hm C^o-|T(o)rt^/-'--r(6)|a.o^) 

Here [r/2] denotes the first integer smaller than or equal to r/2 and we have 
set T = T2N,M- A graphical representation of the transfer matrices appearing 
in the numerator of Eq. (20) is shown in Fig. 7. This new transfer matrix 
can again be treated with the DMRG algorithm described in section 3 where 
either a t or ii plaquette is added corresponding to a decrease in temperature 
T or an increase in real-time t, respectively. 

To demonstrate the method, results for the longitudinal spin-spin auto- 
correlation function of the XXZ-chahi at infinite temperature are shown in 
Fig. 8. For A = the XXZ-model corresponds to free spinless fermions and is 
exactly solvable. We focus on the case of free fermions, as here the analysis of 
the dynamical TMRG (DTMRG) method, its results and numerical errors can 
be done to much greater extent than in the general case. The performance of 
the DTMRG itself is expected to be independent of the strength of the interac- 
tion. The comparison with the exact result in Fig. 8 shows that the maximum 
time before the DTMRG algorithm breaks down increases with the number of 
states. However, the improvement when taking N = 400 instead of = 300 
states is marginal. The reason for the breakdown of the DTMRG computation 
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T(0) T ^''^^-^ T(0) 

Fig. 7. Transfer matrices appearing in the numerator of Eq. (20) for r > 1 with r 
even. The 2 black dots denote the operator O. T, T{0) consist of three parts: A part 
representing exp{—PH) (vertically striped plaquettes), another for exp{itH) (stripes 
from lower left to upper right) and a third part describing exp(— it//) (upper left to 
lower right). T,T{0) are split into system (5*) and environment {E). 




Fig. 8. Autocorrelation function for Z\ = and A — 1 (inset) at T = oo where 
A'' = 50 — 400 states have been kept and S — 0.1. The exact result is shown for 
comparison in the case A = 0. 
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Fig. 9. Largest 100 eigenvalues Ai of ps for Zi = and T = oo calculated exactly. 
The inset shows the discarded weig 

can be traced back to an increase of the discarded weight (see inset of Fig. 9). 
Throughout the RG procedure we keep only N of the leading eigenstates of 
the reduced density matrix ps- As long as the discarded states carry a total 
weight less than, say, 10"'^ the results are faithful. For infinite temperature 
and Zi = we could explain the rapid increase of the discarded weight with 
time by deriving an explicit expression for the leading eigcnstatc of the QTM 
as well as for the corresponding reduced density matrix. At the free fermion 
point the spectrum of this density matrix is multiplicative. Hence, from the 
one-particle spectrum which is calculated by simple numerics we obtain the 
entire spectrum. As shown in Fig. 9 this spectrum becomes more dense with 
increasing time thus setting a characteristic time scale tc{N), quite indepen- 
dent of the discretization 5 of the real time, where the algorithm breaks down. 
Despite these limitations, it is often possible to extrapolate the numerical data 
to larger times using physical arguments thus allowing to obtain frequency- 
dependent quantities by a direct Fourier transform. This way the spin-lattice 
relaxation rate for the Heisenberg chain has been successfully calculated [32] . 
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